Dynamics of woodland soil microbial communities during restoration after industrial and agricultural land use
Bana
2024-06-11
1 Install Packages
# # CRAN
#
# install.packages("vegan")
# install.packages("Polychrome")
# install.packages("dendextend")
# install.packages("ggplotify")
# install.packages("parallelMap")
# install.packages("caret")
# install.packages("randomForest")
#
# # Bioconductor
#
# install.packages("BiocManager")
# BiocManager::install("phyloseq")
# BiocManager::install("ANCOMBC")
# BiocManager::install("mixOmics")
#
# # GitHub
#
# install.packages("devtools")
# devtools::install_github("jbisanz/qiime2R")library(dplyr)
library(ggplot2)2 Load data to Phyloseq
Data is ASVs
# Required library
library(qiime2R)
library(phyloseq)
library(vegan)
# Source files:
feature_table_qza <- "feature_table.qza"
rooted_tree_qza <- "rooted_tree.qza"
taxonomy_qza <- "taxonomy.qza"
metadata_tsv <- "samples.txt"
# Read data
data.phy <- qza_to_phyloseq(
features = feature_table_qza,
tree = rooted_tree_qza,
taxonomy = taxonomy_qza,
metadata = metadata_tsv
)# Cleanup
# Removal of not needed objects, packages and cleaning the RAM
rm(feature_table_qza, rooted_tree_qza, taxonomy_qza, metadata_tsv) # remove unnecessary objects
detach("package:qiime2R", unload=TRUE) # unload qiime2R package
gc() # free unused R memory## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 5098515 272.3 8456080 451.7 NA 8456080 451.7
## Vcells 29988918 228.8 88611965 676.1 16384 91949714 701.63 Filtering and merging
3.1 Filter Industrial and Agriculture samples
# Keep only samples with Industrial or Agricultural prior use
table(sample_data(data.phy)$Former_landuse)##
## Agriculture Ancient_woodland Industrial Rewilding
## 150 20 150 10IndAgri <- subset_samples(data.phy,
Former_landuse %in% c("Industrial","Agriculture"))
IndAgri ## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 60355 taxa and 300 samples ]
## sample_data() Sample Data: [ 300 samples by 11 sample variables ]
## tax_table() Taxonomy Table: [ 60355 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 60355 tips and 60236 internal nodes ]3.2 Filter bacterial Kingdom
table(as.data.frame(tax_table(data.phy))$Kingdom)##
## d__Archaea d__Bacteria Unassigned
## 161 60173 21# Keep only Bacteria
BacIndAgri <- subset_taxa(IndAgri, Kingdom == "d__Bacteria")
BacIndAgri # The OTU count dropped from 60,355 to 60,173## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 60173 taxa and 300 samples ]
## sample_data() Sample Data: [ 300 samples by 11 sample variables ]
## tax_table() Taxonomy Table: [ 60173 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 60173 tips and 60054 internal nodes ]3.3 Merging and Rarefaction
plot_rarefaction_curve <- function(phy.obj, taxa_level) {
# sample depth Curve before Rarefaction
# Extract the ASV table from the phyloseq object
otu_table <- otu_table(phy.obj)
# Transpose the table
if (taxa_are_rows(otu_table)) {
otu_table <- t(otu_table)
}
# Convert to a matrix
otu_matrix <- as(otu_table, "matrix")
# Generate rarefaction curve to decide sample depth
rare_curve <- rarecurve(otu_matrix, step = 100, cex = 0.5, col = "blue", label = FALSE,
xlab = "Sample Size", ylab = taxa_level)
}3.4 1. Family level merging
# Agglomerate Bacteria to family level
BacIndAgriFamily <- tax_glom(BacIndAgri, taxrank="Family")
BacIndAgriFamily## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 769 taxa and 300 samples ]
## sample_data() Sample Data: [ 300 samples by 11 sample variables ]
## tax_table() Taxonomy Table: [ 769 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 769 tips and 768 internal nodes ]3.4.1 Normalization (Rarefaction)
# Check the sample depth for deciding the rarefaction threshold
plot_rarefaction_curve(BacIndAgriFamily, "ASVs (Family)")
3.5 2. Genera level merging
# Agglomerate Bacteria to genus level
BacIndAgriGenera <- tax_glom(BacIndAgri, taxrank="Genus")
BacIndAgriGenera## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1465 taxa and 300 samples ]
## sample_data() Sample Data: [ 300 samples by 11 sample variables ]
## tax_table() Taxonomy Table: [ 1465 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1465 tips and 1464 internal nodes ]3.5.1 Normalization (Rarefaction)
# Check the sample depth for deciding the rarefaction threshold
plot_rarefaction_curve(BacIndAgriGenera, "ASVs (Genus)")
3.6 3. Species level merging
3.6.1 Normalization (Rarefaction)
# Check the sample depth for deciding the rarefaction threshold
plot_rarefaction_curve(BacIndAgri, "ASVs (Species)")
3.7 Applying Rarefaction on Genera merged data
# Rarefy the phyloseq object to an even depth of 10000 sequences per sample
BacIndAgriGenera.rarefied <- rarefy_even_depth(BacIndAgriGenera,
sample.size = 10000,
rngseed = 123, replace = FALSE)Rarefation done on Genera level, for multiple reasons: 1. not all ASVs are classified up to species 2. The 16srRNA sequencing achieves better genus-level resolution
The decision of sample depth threshold (10000) to use for normalization was mainly based on visual assessment of two things: sample reaching the plateau of ASVs count and avoiding losing any samples.
3.8 Filtering low quality ASVs (at Genera level)
# (Function) to create and save a histogram
save_histogram <- function(data, xlab, ylab, main, filename) {
# Display the histogram
hist(data, xlab = xlab, ylab = ylab, main = main)
# Save the histogram to a file
dev.copy(png, filename)
dev.off()
}# Convert the OTU table to matrix
BacIndAgriGenera.numeric <- as.numeric(otu_table(BacIndAgriGenera))
BacIndAgriGenera.mtx <- matrix(BacIndAgriGenera.numeric, nrow=nrow(otu_table(BacIndAgriGenera)))
# Getting the ASV that has more than 0 count in each sample THRESHOLD
BacIndAgriGenera.count_filter <- BacIndAgriGenera.mtx > 0
# Sum the number of ASV for each sample that has > count threshold
# The result will show the number of ASV (frequency, y-axis) > count threshold in each sample (x-axis)
BacIndAgriGenera.count_filter.sum <- apply(BacIndAgriGenera.count_filter, 1, sum)
# Display and save the first histogram
save_histogram(BacIndAgriGenera.count_filter.sum, xlab = "Samples", ylab = "ASVs count (Genera)",
main = "Frequency of ASVs with count > 0 across samples",
filename = "results/filter_low_quality_ASVs/ASVsGenera_count_histogram.png")
## quartz_off_screen
## 2# Plot histogram for sample interval (0-10)
BacIndAgriGenera.count_filter.sum.1_10.interval <- BacIndAgriGenera.count_filter.sum[BacIndAgriGenera.count_filter.sum < 10]
# Display and save the second histogram
save_histogram(BacIndAgriGenera.count_filter.sum.1_10.interval, xlab = "Samples", ylab = "ASVs count (Genera)",
main = "Frequency of ASVs with count > 0 in 10 samples or less",
filename = "results/filter_low_quality_ASVs/ASVsGenera_count_histogram_1_10_interval.png")
## quartz_off_screen
## 2Interpretation: Let’s say for the second plot we have ~ 150 ASVs that are present in two samples.
all_taxa.num <- length(BacIndAgriGenera.count_filter.sum)
taxa_less10samples <- length(BacIndAgriGenera.count_filter.sum.1_10.interval)
lost_taxa_perc <- round((taxa_less10samples / all_taxa.num) * 100, 2)
print(paste("Out of", all_taxa.num, "the number of lost ASVs (Genera) after filtering taxa with abundance > 0 counts in less than 10 samples", taxa_less10samples, "Which is", lost_taxa_perc, "%"))## [1] "Out of 1465 the number of lost ASVs (Genera) after filtering taxa with abundance > 0 counts in less than 10 samples 766 Which is 52.29 %"From count filtering histogram and stats, half of ASVs (Genera) will be lost. Therefore, for further analysis the data will be kept without filtering unless it’s required for a statistical method to be used.
4 Visualizing Microbial diversity
library(ggplot2)
library(plotly)
library(RColorBrewer)4.1 Bar plot
4.1.1 Former land use
# Aggregate data at Phylum level
BacIndAgriPhylum.rarefied <- tax_glom(BacIndAgriGenera.rarefied, taxrank = "Phylum")4.2 Check data types
# Make sure pH / Conductivity / Age are numerical
sample_data(BacIndAgriPhylum.rarefied)$pH <- as.numeric(as.character(sample_data(BacIndAgriPhylum.rarefied)$pH))
sample_data(BacIndAgriPhylum.rarefied)$Conductivity <- as.numeric(as.character(sample_data(BacIndAgriPhylum.rarefied)$Conductivity))
sample_data(BacIndAgriPhylum.rarefied)$Woodland_age <- as.numeric(as.character(sample_data(BacIndAgriPhylum.rarefied)$Woodland_age))env_factor_bar_plot <- function(phyloseq.obj, env_factor, save_path, levels=c()) {
# Merge samples by environmental factor (categorical one)
phyloseq.obj.merged <- merge_samples(phyloseq.obj, env_factor)
# Transform counts to relative abundances
phyloseq.obj.merged <- transform_sample_counts(phyloseq.obj.merged, function(x) x / sum(x))
# Extract the data for plotting
phyloseq.obj.merged.df <- psmelt(phyloseq.obj.merged)
# Get counts per phylum
Phyla.df <- phyloseq.obj.merged.df %>%
group_by(Phylum) %>%
summarise(Count = sum(Abundance))
# Select the cut-off for the Phylum taxa (e.g. 1% of total count)
cutoff <- 0.01 * sum(phyloseq.obj.merged.df$Abundance)
# Select low-abundant Phyla (with total counts below the cutoff)
lowAbundant <- Phyla.df[Phyla.df$Count <= cutoff,]$Phylum
# Substitute Phylum names to "<1%" for the low-abundant phyla
phyloseq.obj.merged.df[phyloseq.obj.merged.df$Phylum %in% lowAbundant,]$Phylum <- '<1%'
# Ensure env_factor levels are ordered correctly in the melted data frame
if (length(levels) > 0){
phyloseq.obj.merged.df$Sample <- factor(phyloseq.obj.merged.df$Sample, levels = levels)
}
# Define high contrast palette for bar plot
n_colors <- length(unique(phyloseq.obj.merged.df$Phylum))
high_contrast_palette <- c(brewer.pal(min(9, n_colors), "Set1"), brewer.pal(max(0, n_colors - 9), "Dark2"))
# Create bar plot
barplot <- ggplot(phyloseq.obj.merged.df, aes(x = Sample, y = Abundance, fill = Phylum, text = paste("Phylum:", Phylum, "<br>", env_factor, ":", Sample, "<br>Relative Abundance:", Abundance))) +
geom_bar(stat = "identity", position = "stack") +
scale_fill_manual(values = high_contrast_palette) +
theme(axis.text.x = element_text(angle = 0, hjust = 1)) +
labs(title = paste("Taxa Abundance by", env_factor), x = env_factor, y = "Relative Abundance")
# Convert to plotly
barplot.plotly <- ggplotly(barplot, tooltip = "text")
barplot.plotly
# Save the plot as an HTML file
htmlwidgets::saveWidget(barplot.plotly, save_path)
return(barplot.plotly)
} barplot.FLU <- env_factor_bar_plot(phyloseq.obj = BacIndAgriPhylum.rarefied, env_factor = "Former_landuse",
save_path = "results/plotly/Taxonomy_Phyla_landuse_normalized.html")
barplot.FLU4.2.1 pH
pHValues <- as.numeric(sample_data(BacIndAgriGenera)$pH)## Warning: NAs introduced by coercionpHValues.removed_nan <- na.omit(pHValues)
print("")## [1] ""print("max pH")## [1] "max pH"max(pHValues.removed_nan)## [1] 7.98print("Min pH")## [1] "Min pH"min(pHValues.removed_nan)## [1] 3.78print("Median pH")## [1] "Median pH"median(pHValues.removed_nan)## [1] 5.54# Make sure pH column is numeric
sample_data(BacIndAgriPhylum.rarefied)$pH <- as.numeric(as.character(sample_data(BacIndAgriPhylum.rarefied)$pH))
# Categorizing pH
sample_data(BacIndAgriPhylum.rarefied)$pH_category <- factor(
cut(sample_data(BacIndAgriPhylum.rarefied)$pH,
breaks = c(3.78, 5, 6, 7.98),
labels = c("<5", "5-6", ">6")),
levels = c("<5", "5-6", ">6")
)pH_less_5 <- as.numeric(sample_data(BacIndAgriPhylum.rarefied)$pH_category == "<5")
pH_less_5 <- na.omit(pH_less_5)
print(paste("Number of samples with to pH < 5 =", sum(pH_less_5)))## [1] "Number of samples with to pH < 5 = 50"pH_5_6 <- as.numeric(sample_data(BacIndAgriPhylum.rarefied)$pH_category == "5-6")
pH_5_6 <- na.omit(pH_5_6)
print(paste("Number of samples with to pH 5-6 =", sum(pH_5_6)))## [1] "Number of samples with to pH 5-6 = 145"pH_more_6 <- as.numeric(sample_data(BacIndAgriPhylum.rarefied)$pH_category == ">6")
pH_more_6 <- na.omit(pH_more_6)
print(paste("Number of samples with pH >6 =", sum(pH_more_6)))## [1] "Number of samples with pH >6 = 95"barplot.pH <- env_factor_bar_plot(phyloseq.obj = BacIndAgriPhylum.rarefied, env_factor = "pH_category",
save_path = "results/plotly/Taxonomy_Phyla_pH_normalized.html",
levels = c("<5", "5-6", ">6"))
barplot.pH4.2.2 Conductivity
conductivityValues <- as.numeric(sample_data(BacIndAgriGenera)$Conductivity)## Warning: NAs introduced by coercion# Calculate the quartiles
quartiles <- quantile(conductivityValues, probs = c(0.25, 0.5, 0.75), na.rm = TRUE)
quartiles## 25% 50% 75%
## 39 64 134# Make sure conductivity column is numeric
sample_data(BacIndAgriPhylum.rarefied)$Conductivity <- as.numeric(as.character(sample_data(BacIndAgriPhylum.rarefied)$Conductivity))
# Calculate the quartiles for conductivity
quartiles <- quantile(sample_data(BacIndAgriPhylum.rarefied)$Conductivity, probs = c(0.25, 0.5, 0.75), na.rm = TRUE)
# Categories based on quartiles
sample_data(BacIndAgriPhylum.rarefied)$Conductivity_category <-
cut(sample_data(BacIndAgriPhylum.rarefied)$Conductivity,
breaks = c(-Inf, quartiles, Inf),
labels = c("Q1: <0.25", "Q2: 0.25-0.5", "Q3: 0.5-0.75", "Q4: 0.75<"),
include.lowest = TRUE)barplot.conductivity <- env_factor_bar_plot(phyloseq.obj = BacIndAgriPhylum.rarefied, env_factor = "Conductivity_category",
save_path = "results/plotly/Taxonomy_Phyla_Conductivity_normalized.html")
barplot.conductivity4.3 Line Plot
4.3.1 FLU vs pH
# Make sure pH column is numeric
sample_data(BacIndAgri)$pH <- as.numeric(as.character(sample_data(BacIndAgri)$pH))
# Make Former_landuse column as factor
sample_data(BacIndAgri)$Former_landuse <- as.factor(sample_data(BacIndAgri)$Former_landuse)
# Extract sample data into a df
FLU_pH.df <- data.frame(sample_data(BacIndAgri))
FLU_pH.plot <- plot_ly()
# Compute density estimates for each Former_landuse level
for(level in levels(FLU_pH.df$Former_landuse)) {
data_subset <- FLU_pH.df[FLU_pH.df$Former_landuse == level, ]
dens <- density(data_subset$pH, na.rm = TRUE) # Compute density
# Add density trace
FLU_pH.plot <- FLU_pH.plot %>% add_trace(
x = dens$x,
y = dens$y,
type = 'scatter',
mode = 'lines',
name = level,
opacity = 0.6
)
}
# Customize the layout
FLU_pH.plot <- FLU_pH.plot %>% layout(
title = "Kernel Density Estimate of pH by Former Landuse",
xaxis = list(title = "pH"),
yaxis = list(title = "Density")
)
FLU_pH.plot# Save the plot as an HTML file
htmlwidgets::saveWidget(FLU_pH.plot, "results/plotly/pH_vs_FLU.html")4.4 Correlation
Scatter Plot ### pH vs Conductivity
library(broom)# Make sure conductivity column is numeric
sample_data(BacIndAgri)$Conductivity <- as.numeric(as.character(sample_data(BacIndAgri)$Conductivity))## Warning: NAs introduced by coercion# Extract sample data
Cond_pH.df <- data.frame(sample_data(BacIndAgri))
# linear regression
model <- lm(Conductivity ~ pH, data = Cond_pH.df)
# Extract the coefficients
summary_lm <- summary(model) # Benjamini-Hochberg adjusted p-value by default
print(summary_lm)##
## Call:
## lm(formula = Conductivity ~ pH, data = Cond_pH.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -85.67 -37.37 -11.57 29.63 237.12
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -96.885 20.258 -4.783 2.74e-06 ***
## pH 31.986 3.474 9.207 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 56.62 on 293 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.2244, Adjusted R-squared: 0.2218
## F-statistic: 84.77 on 1 and 293 DF, p-value: < 2.2e-16# Extract R-squared and slope
r_squared <- summary_lm$r.squared
slope <- summary_lm$coefficients[2, 1]
# Create the regression plot
Cond_pH.plot <- ggplot(Cond_pH.df, aes(x = pH, y = Conductivity, color = Former_landuse)) +
geom_point() +
geom_smooth(method = "lm", se = TRUE) +
labs(title = "Regression of Conductivity on pH by Former Land Use",
subtitle = paste("R-squared =", round(r_squared, 2),
"Slope =", round(slope, 2)),
x = "pH",
y = "Conductivity") +
theme_minimal()
print(Cond_pH.plot)## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
# Convert to plotly
Cond_pH.plot_plotly <- ggplotly(Cond_pH.plot)## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_smooth()`).print(Cond_pH.plot_plotly)
# Save the plot as an HTML file
htmlwidgets::saveWidget(Cond_pH.plot_plotly, "results/plotly/Regression_pH_Conductivity_by_LandUse.html")Adjusted R-squared 0.2218 (very low) shows that the model failed to explain the real data points and the results are not reliable. For the results it tells that for one unit increase of pH the conductivity increases by 31.986, and theoretically at pH 0 the conductivity is -96.885
4.4.1 Woodland age vs taxa abundance
library(dplyr)# Function to compute the correlation between microbial abundance and environmental factor
abundance_env_factor_corr <- function(phyloseq.obj, env_factor){
# Create abundance df
phyloseq.obj.abund <- psmelt(phyloseq.obj)
# Ensure environmental factors are numerical
phyloseq.obj.abund$Woodland_age <- as.numeric(as.character(phyloseq.obj.abund$Woodland_age))
phyloseq.obj.abund$Conductivity <- as.numeric(as.character(phyloseq.obj.abund$Conductivity))
phyloseq.obj.abund$pH <- as.numeric(as.character(phyloseq.obj.abund$pH))
# Aggregate abundance by Site_code (as samples of each site share the same value of pH, Conductivity, age) and Phylum, summing ASVs
phyloseq.obj.abund.summary <- phyloseq.obj.abund %>%
group_by(Site_code, Phylum) %>%
summarise(Abundance = mean(Abundance), .groups = 'drop')
# Align Sample IDs between phyloseq sample data and the summary data
sample_data_df <- as.data.frame(sample_data(phyloseq.obj))
sample_data_df$Sample <- rownames(sample_data_df)
# Join the summary data with the sample metadata
phyloseq.obj.abund.summary <- left_join(phyloseq.obj.abund.summary, sample_data_df, by = "Site_code")
# Correlation function
compute_correlation <- function(df) {
cor_result <- cor.test(df[[env_factor]], df$Abundance, method = "spearman")
data.frame(
corr = cor_result$estimate,
p_value = cor_result$p.value
)
}
# Compute correlation and p-value for each Phylum
correlations <- phyloseq.obj.abund.summary %>%
group_by(Phylum) %>%
summarise(correlation = list(compute_correlation(cur_data())), .groups = 'drop')
# Flatten the list column to extract correlation and p-value
correlations <- correlations %>%
mutate(
correlation = sapply(correlations$correlation, function(x) x$corr),
p_value = sapply(correlations$correlation, function(x) x$p_value)
)
# Calculate the adjusted p-value
correlations$pAdjust <- p.adjust(correlations$p_value)
# Filter taxa based on adjusted p-value <= 0.05
significant_taxa <- correlations %>%
filter(pAdjust <= 0.05)
# Get top positively and negatively correlated taxa
top_positive <- significant_taxa %>%
arrange(desc(correlation)) %>%
slice_head(n = 4)
top_negative <- significant_taxa %>%
arrange(correlation) %>%
slice_head(n = 4)
print("Top positively correlated taxa with Woodland age:")
print(top_positive)
print("Top negatively correlated taxa with Woodland age:")
print(top_negative)
return(list("abund.summary"=phyloseq.obj.abund.summary,
"top_positive"=top_positive, "top_negative"=top_negative))
}age_abund_corr <- abundance_env_factor_corr(phyloseq.obj = BacIndAgriPhylum.rarefied,
env_factor = "Woodland_age")## Warning in left_join(phyloseq.obj.abund.summary, sample_data_df, by = "Site_code"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 146 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Warning: There were 58 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `correlation = list(compute_correlation(cur_data()))`.
## ℹ In group 1: `Phylum = "Acidobacteriota"`.
## Caused by warning:
## ! `cur_data()` was deprecated in dplyr 1.1.0.
## ℹ Please use `pick()` instead.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 57 remaining warnings.## [1] "Top positively correlated taxa with Woodland age:"
## # A tibble: 4 × 4
## Phylum correlation p_value pAdjust
## <chr> <dbl> <dbl> <dbl>
## 1 CSP1-3 0.370 3.70e-11 0.00000000211
## 2 Sumerlaeota 0.304 7.64e- 8 0.00000428
## 3 Desulfobacterota_C 0.299 1.29e- 7 0.00000707
## 4 KSB1 0.283 6.44e- 7 0.0000348
## [1] "Top negatively correlated taxa with Woodland age:"
## # A tibble: 4 × 4
## Phylum correlation p_value pAdjust
## <chr> <dbl> <dbl> <dbl>
## 1 Firmicutes_E -0.278 0.00000104 0.0000552
## 2 Firmicutes_A -0.261 0.00000468 0.000234
## 3 Firmicutes_B_370531 -0.222 0.000106 0.00454
## 4 Planctomycetota -0.199 0.000520 0.0218top_pos_neg_plot <- function(top_positive, top_negative, abund.summary, env_factor) {
# Select the top positively and negatively correlated taxa
top_taxa <- rbind(top_positive[1, ], top_negative[1, ])
# Extract data for these taxa
top_taxa.abund.summary <- abund.summary %>%
filter(Phylum %in% top_taxa$Phylum)
# Scatter plots
for (phylum in unique(top_taxa$Phylum)) {
# Filter data for the current Phylum
phylum_data <- top_taxa.abund.summary %>% filter(Phylum == phylum)
# Fit a linear model
lm_formula <- paste("Abundance", "~", env_factor)
lm_formula <- as.formula(lm_formula)
lm_model <- lm(lm_formula, data = phylum_data)
summary_lm <- summary(lm_model)
print("--------------------------------------")
print(summary_lm)
# Extract R-squared and slope
r_squared <- summary_lm$r.squared
slope <- summary_lm$coefficients[2, 1]
# Create scatter plot
p <- ggplot(phylum_data, aes_string(x = env_factor, y = "Abundance")) +
geom_point(color = "blue") +
geom_smooth(method = "lm", se = TRUE, color = "red") +
labs(title = paste("Abundance of", phylum, "vs", env_factor),
subtitle = paste("R-squared =", round(r_squared, 2),
"Slope =", round(slope, 2)),
x = env_factor,
y = "Abundance") +
theme_minimal()
# Print the plot
print(p)
}
}top_pos_neg_plot(top_positive = age_abund_corr$top_positive, top_negative = age_abund_corr$top_negative,
abund.summary = age_abund_corr$abund.summary, env_factor = "Woodland_age")## [1] "--------------------------------------"
##
## Call:
## lm(formula = lm_formula, data = phylum_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.27042 -0.10884 -0.02256 0.03078 0.99233
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.150011 0.034941 -4.293 2.38e-05 ***
## Woodland_age 0.006275 0.000851 7.373 1.65e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2391 on 298 degrees of freedom
## Multiple R-squared: 0.1543, Adjusted R-squared: 0.1515
## F-statistic: 54.37 on 1 and 298 DF, p-value: 1.654e-12## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `geom_smooth()` using formula = 'y ~ x'## [1] "--------------------------------------"
##
## Call:
## lm(formula = lm_formula, data = phylum_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5679 -0.3770 -0.2071 -0.0433 4.7929
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.696519 0.117930 5.906 9.54e-09 ***
## Woodland_age -0.008039 0.002872 -2.799 0.00547 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.807 on 298 degrees of freedom
## Multiple R-squared: 0.02561, Adjusted R-squared: 0.02234
## F-statistic: 7.832 on 1 and 298 DF, p-value: 0.005467## `geom_smooth()` using formula = 'y ~ x'
4.4.2 pH vs taxa abundance
pH_abund_corr <- abundance_env_factor_corr(phyloseq.obj = BacIndAgriPhylum.rarefied,
env_factor = "pH")## Warning in left_join(phyloseq.obj.abund.summary, sample_data_df, by = "Site_code"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 146 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Warning: There were 57 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `correlation = list(compute_correlation(cur_data()))`.
## ℹ In group 1: `Phylum = "Acidobacteriota"`.
## Caused by warning in `cor.test.default()`:
## ! Cannot compute exact p-value with ties
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 56 remaining warnings.## [1] "Top positively correlated taxa with Woodland age:"
## # A tibble: 4 × 4
## Phylum correlation p_value pAdjust
## <chr> <dbl> <dbl> <dbl>
## 1 Chloroflexota 0.699 1.34e-44 7.65e-43
## 2 Methylomirabilota 0.506 1.53e-20 8.10e-19
## 3 Myxococcota_A_437813 0.438 3.02e-15 1.51e-13
## 4 Myxococcota_A_473307 0.425 2.25e-14 1.10e-12
## [1] "Top negatively correlated taxa with Woodland age:"
## # A tibble: 4 × 4
## Phylum correlation p_value pAdjust
## <chr> <dbl> <dbl> <dbl>
## 1 Eremiobacterota -0.678 5.35e-41 2.99e-39
## 2 Planctomycetota -0.675 1.23e-40 6.78e-39
## 3 Acidobacteriota -0.559 1.17e-25 6.31e-24
## 4 Firmicutes_D -0.505 1.58e-20 8.20e-19top_pos_neg_plot(top_positive = pH_abund_corr$top_positive, top_negative = pH_abund_corr$top_negative,
abund.summary = pH_abund_corr$abund.summary, env_factor = "pH")## [1] "--------------------------------------"
##
## Call:
## lm(formula = lm_formula, data = phylum_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -273.56 -95.68 -8.36 69.29 497.47
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -200.772 54.059 -3.714 0.000244 ***
## pH 106.065 9.271 11.441 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 151.1 on 293 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.3088, Adjusted R-squared: 0.3064
## F-statistic: 130.9 on 1 and 293 DF, p-value: < 2.2e-16## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).## [1] "--------------------------------------"
##
## Call:
## lm(formula = lm_formula, data = phylum_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.705 -8.128 -2.338 5.547 48.882
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 58.6198 4.1269 14.20 <2e-16 ***
## pH -8.2371 0.7077 -11.64 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.53 on 293 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.3162, Adjusted R-squared: 0.3138
## F-statistic: 135.5 on 1 and 293 DF, p-value: < 2.2e-16## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 5 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
4.4.3 Conductivity vs taxa abundance
conductivity_abund_corr <- abundance_env_factor_corr(phyloseq.obj = BacIndAgriPhylum.rarefied,
env_factor = "Conductivity")## Warning in left_join(phyloseq.obj.abund.summary, sample_data_df, by = "Site_code"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 146 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Warning: There were 57 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `correlation = list(compute_correlation(cur_data()))`.
## ℹ In group 1: `Phylum = "Acidobacteriota"`.
## Caused by warning in `cor.test.default()`:
## ! Cannot compute exact p-value with ties
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 56 remaining warnings.## [1] "Top positively correlated taxa with Woodland age:"
## # A tibble: 4 × 4
## Phylum correlation p_value pAdjust
## <chr> <dbl> <dbl> <dbl>
## 1 Actinobacteriota 0.699 1.69e-44 9.66e-43
## 2 Chloroflexota 0.546 2.66e-24 1.33e-22
## 3 Firmicutes_B_370539 0.454 2.16e-16 9.94e-15
## 4 Methylomirabilota 0.415 9.78e-14 4.40e-12
## [1] "Top negatively correlated taxa with Woodland age:"
## # A tibble: 4 × 4
## Phylum correlation p_value pAdjust
## <chr> <dbl> <dbl> <dbl>
## 1 Acidobacteriota -0.678 4.25e-41 2.38e-39
## 2 FCPU426 -0.677 5.91e-41 3.25e-39
## 3 Elusimicrobiota -0.639 3.03e-35 1.64e-33
## 4 Armatimonadota -0.632 2.47e-34 1.31e-32top_pos_neg_plot(top_positive = conductivity_abund_corr$top_positive,
top_negative = conductivity_abund_corr$top_negative,
abund.summary = conductivity_abund_corr$abund.summary, env_factor = "Conductivity")## [1] "--------------------------------------"
##
## Call:
## lm(formula = lm_formula, data = phylum_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2040.42 -602.79 -96.51 430.54 2617.57
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1432.4810 80.4577 17.80 <2e-16 ***
## Conductivity 9.3798 0.7438 12.61 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 818.5 on 293 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.3518, Adjusted R-squared: 0.3496
## F-statistic: 159 on 1 and 293 DF, p-value: < 2.2e-16## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).## [1] "--------------------------------------"
##
## Call:
## lm(formula = lm_formula, data = phylum_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1269.72 -351.81 -18.32 318.38 1465.93
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2601.1759 52.4365 49.61 <2e-16 ***
## Conductivity -5.5302 0.4848 -11.41 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 533.5 on 293 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.3075, Adjusted R-squared: 0.3052
## F-statistic: 130.1 on 1 and 293 DF, p-value: < 2.2e-16## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 5 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
5 Alpha Diversity
5.1 plot
alpha_diversity.plot <- function(alpha_richness, phyloseq_obj, env_factor){
# Apply Wilcoxon test
chao1_res <- wilcox.test(alpha_richness$Chao1~
sample_data(phyloseq_obj)[[env_factor]])
chao1.p_value <- chao1_res$p.value
shannon_res <- wilcox.test(alpha_richness$Shannon~
sample_data(phyloseq_obj)[[env_factor]])
shannon.p_value <- shannon_res$p.value
simpson_res <- wilcox.test(alpha_richness$Simpson~
sample_data(phyloseq_obj)[[env_factor]])
simpson.p_value <- simpson_res$p.value
plot_richness(phyloseq_obj,
measures=c("Chao1", "Shannon", "Simpson"),
color=env_factor, x=env_factor) +
geom_boxplot() + ggtitle(paste("Alpha diversity vs", env_factor),
subtitle=paste("Chao1 p-value:", chao1.p_value,
"\nShannon p-value:", shannon.p_value,
"\nSimpson p-value:", simpson.p_value))
}5.1.1 by FLU
# Calculate Alpha diversity indices
alpha_richness = estimate_richness(
BacIndAgriGenera.rarefied, measures = c("Chao1", "Shannon", "Simpson"))alpha_diversity.plot(alpha_richness = alpha_richness, phyloseq_obj = BacIndAgriGenera.rarefied,
env_factor = "Former_landuse")
### by Region
alpha_diversity.plot(alpha_richness = alpha_richness, phyloseq_obj = BacIndAgriGenera.rarefied,
env_factor = "Region")
5.2 Variance
5.2.1 Shannon
shannon_values <- alpha_richness$Shannon
site_codes <- sample_data(BacIndAgriGenera.rarefied)$Site_code
FLU_values <- sample_data(BacIndAgriGenera.rarefied)$Former_landuse
shannon_site.df <- data.frame(Site_code = site_codes, Shannon = shannon_values,
Former_landuse = FLU_values)ggplot(shannon_site.df, aes(x = Site_code, y = Shannon, fill = Former_landuse)) +
geom_boxplot() +
theme_minimal() +
labs(title = "Shannon Diversity by Site",
x = "Site",
y = "Shannon Index") +
theme(axis.text.x = element_text(angle = 90, size=4)) 
# Calculate variance for sites replicas
variance_data.shannon <- shannon_site.df %>%
group_by(Site_code) %>%
summarise(CV = (sd(Shannon, na.rm = TRUE)/mean(Shannon, na.rm = TRUE))*100, Former_landuse = unique(Former_landuse))ggplot(variance_data.shannon, aes(x = Site_code, y = CV, fill = Former_landuse)) +
geom_bar(stat = "identity") +
theme_minimal() +
labs(title = "Coefficient of Variation of Shannon Diversity by Site",
x = "Site",
y = "Coefficient of Variation") +
theme(axis.text.x = element_text(angle = 90, size=4)) 
library(readxl)sites_pH_EC <- read_excel("pH_EC_RestREco_8th_July_2024.xlsx")5.2.2 pH
ggplot(sites_pH_EC, aes(x = Site_code, y = pH, fill = Former_landuse)) +
geom_boxplot() +
theme_minimal() +
labs(title = "pH Diversity by Site",
x = "Site",
y = "pH") +
theme(axis.text.x = element_text(angle = 90, size = 6)) 
# Calculate variance for sites replicas
variance_data.pH <- sites_pH_EC %>%
group_by(Site_code) %>%
summarise(CV = (sd(pH, na.rm = TRUE)/mean(pH, na.rm = TRUE))*100, Former_landuse = unique(Former_landuse))ggplot(variance_data.pH, aes(x = Site_code, y = CV, fill = Former_landuse)) +
geom_bar(stat = "identity") +
theme_minimal() +
labs(title = "Coefficient of Variation of pH by Site",
x = "Site",
y = "Coefficient of Variation") +
theme(axis.text.x = element_text(angle = 90, size = 6)) 
5.2.3 Conductivity
ggplot(sites_pH_EC, aes(x = Site_code, y = Conductivity, fill = Former_landuse)) +
geom_boxplot() +
theme_minimal() +
labs(title = "Conductivity Diversity by Site",
x = "Site",
y = "Conductivity") +
theme(axis.text.x = element_text(angle = 90, size = 6)) 
# Calculate variance for sites replicas
variance_data.conductivity <- sites_pH_EC %>%
group_by(Site_code) %>%
summarise(CV = (sd(Conductivity, na.rm = TRUE)/mean(Conductivity, na.rm = TRUE))*100,
Former_landuse = unique(Former_landuse))ggplot(variance_data.conductivity, aes(x = Site_code, y = CV, fill = Former_landuse)) +
geom_bar(stat = "identity") +
theme_minimal() +
labs(title = "Coefficient of Variation of Conductivity by Site",
x = "Site",
y = "Coefficient of Variation") +
theme(axis.text.x = element_text(angle = 90, size = 6)) 
5.3 Modelling (nested design)
#install.packages("remotes")
#remotes::install_github("mikemc/speedyseq")library(speedyseq)# Merge replicas for each site by median (non-parametric)
BacIndAgriGenera.rarefied.merged_replica <- merge_samples2(x = BacIndAgriGenera.rarefied, group = "Site_code", fun_otu = median)# Sanity check of counts per sample (the sum should be same across all samples)
plot(colSums(otu_table(BacIndAgriGenera.rarefied.merged_replica)),
ylab="Count", xlab="Samples", xaxt="n",
main="Counts per sample after merging")
# Plot the ASVs count per site after median merging to decide a threshold for normalization
plot_rarefaction_curve(BacIndAgriGenera.rarefied.merged_replica, "ASVs (Genera)")
# Rarefy (normalize) to an even depth of 8000 sequences per sample
BacIndAgriGenera.rarefied2.merged_replica <- rarefy_even_depth(BacIndAgriGenera.rarefied.merged_replica,
sample.size = 8000,
rngseed = 123, replace = FALSE)# Sanity check
plot(colSums(otu_table(BacIndAgriGenera.rarefied2.merged_replica)),
ylab="Count", xlab="Samples", xaxt="n",
main="Counts per sample after merging and rarefaction")
# Calculate Shannon
richness.shannon = estimate_richness(
BacIndAgriGenera.rarefied2.merged_replica, measures = c("Shannon"))
# Add shannon to the phyloseq obj
sample_data(BacIndAgriGenera.rarefied2.merged_replica)$Shannon <- richness.shannon$Shannon# Get sample data into dataframe
BacIndAgriGenera.rarefied2.merged_replica.df <-
as(sample_data(BacIndAgriGenera.rarefied2.merged_replica), "data.frame")reg_model.alphaDiv <- lm(Shannon~Region/Former_landuse, data=BacIndAgriGenera.rarefied2.merged_replica.df)
summary(reg_model.alphaDiv)##
## Call:
## lm(formula = Shannon ~ Region/Former_landuse, data = BacIndAgriGenera.rarefied2.merged_replica.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.33604 -0.11658 0.02269 0.09271 0.40594
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.44821 0.03803 116.951 < 2e-16
## RegionScotland -0.18271 0.05379 -3.397 0.00126
## RegionEngland:Former_landuseIndustrial 0.06363 0.05379 1.183 0.24181
## RegionScotland:Former_landuseIndustrial 0.29169 0.05379 5.423 1.29e-06
##
## (Intercept) ***
## RegionScotland **
## RegionEngland:Former_landuseIndustrial
## RegionScotland:Former_landuseIndustrial ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1473 on 56 degrees of freedom
## Multiple R-squared: 0.3782, Adjusted R-squared: 0.3449
## F-statistic: 11.36 on 3 and 56 DF, p-value: 6.367e-066 Differential Abundance
On merged sites replicas
library(ANCOMBC)Differential Abundance Tests: pH, conductivity, Region, age, FLU
pH + conductivity + age –> Trend test Region –> Global test / pairwise FLU –> pairwise
6.1 FLU
# Estimate differential abundances
# 1. FLU (Pairwise) bonferroni correction is often applied to control for the increased risk of Type I errors
ancom.FLU = ancombc2(data = BacIndAgriGenera.rarefied2.merged_replica, tax_level = "Genus",
fix_formula = "Former_landuse",
group = "Former_landuse",
p_adj_method = "bonferroni", alpha = 0.05, pairwise = TRUE,
prv_cut = 0.10, lib_cut = 0,
n_cl = 4, verbose = TRUE)
# alpha: significance cut-off
# prv_cut: filter taxa not present in prv_cut * samples_number
# lib_cut: filter samples with library size threshold# Extract the results table
ancom.FLU.df <- ancom.FLU$res
colnames(ancom.FLU.df)## [1] "taxon" "lfc_(Intercept)"
## [3] "lfc_Former_landuseIndustrial" "se_(Intercept)"
## [5] "se_Former_landuseIndustrial" "W_(Intercept)"
## [7] "W_Former_landuseIndustrial" "p_(Intercept)"
## [9] "p_Former_landuseIndustrial" "q_(Intercept)"
## [11] "q_Former_landuseIndustrial" "diff_(Intercept)"
## [13] "diff_Former_landuseIndustrial" "passed_ss_(Intercept)"
## [15] "passed_ss_Former_landuseIndustrial"# Select only columns that we need
ancom.FLU.df <- ancom.FLU.df %>%
select(taxon, contains("Former_landuse"))# Select differentially abundant taxa
dif.FLU.df <- ancom.FLU.df %>%
filter(diff_Former_landuseIndustrial & passed_ss_Former_landuseIndustrial) %>%
arrange(desc(lfc_Former_landuseIndustrial))# Filter out values that are not 2 fold change in both positive and negative direction
dif.FLU.df <- dif.FLU.df %>%
filter(lfc_Former_landuseIndustrial >= 1 | lfc_Former_landuseIndustrial <= -1)
# Create the 'direct' column to categorize LFC
dif.FLU.df <- dif.FLU.df %>%
mutate(
direct = ifelse(lfc_Former_landuseIndustrial >= 1, "Positive LFC", "Negative LFC")
)
# Ensure direct factorized
dif.FLU.df$direct =
factor(dif.FLU.df$direct, levels = c("Positive LFC", "Negative LFC"))
dif.FLU.df$taxon =
factor(dif.FLU.df$taxon, levels = dif.FLU.df$taxon)dif.FLU.df## taxon lfc_Former_landuseIndustrial
## 1 Granulicella_B 1.852967
## 2 VAYN01 1.294547
## 3 Sphingomonas_L_486704 1.136363
## 4 Georgfuchsia -1.037880
## 5 Bog-159 -1.056645
## 6 Gp13-AA74 -1.142105
## 7 Clostridium_AD -1.150575
## 8 AV80 -1.333078
## 9 Bacillus_A -1.345840
## 10 Sporosarcina -1.389808
## 11 Bacillus_BK -1.465648
## se_Former_landuseIndustrial W_Former_landuseIndustrial
## 1 0.1366279 13.562139
## 2 0.2358824 5.488104
## 3 0.1619448 7.016978
## 4 0.1380735 -7.516871
## 5 0.2427298 -4.353176
## 6 0.1415131 -8.070667
## 7 0.2120892 -5.424958
## 8 0.1850792 -7.202746
## 9 0.1663562 -8.090111
## 10 0.2340446 -5.938220
## 11 0.1507787 -9.720521
## p_Former_landuseIndustrial q_Former_landuseIndustrial
## 1 9.970687e-06 4.127864e-03
## 2 1.063179e-06 4.401563e-04
## 3 6.082429e-06 2.518126e-03
## 4 5.881795e-07 2.435063e-04
## 5 5.884443e-05 2.436160e-02
## 6 3.237361e-07 1.340267e-04
## 7 1.230925e-06 5.096030e-04
## 8 2.127503e-09 8.807864e-07
## 9 1.914571e-08 7.926325e-06
## 10 1.157035e-06 4.790125e-04
## 11 2.578196e-05 1.067373e-02
## diff_Former_landuseIndustrial passed_ss_Former_landuseIndustrial
## 1 TRUE TRUE
## 2 TRUE TRUE
## 3 TRUE TRUE
## 4 TRUE TRUE
## 5 TRUE TRUE
## 6 TRUE TRUE
## 7 TRUE TRUE
## 8 TRUE TRUE
## 9 TRUE TRUE
## 10 TRUE TRUE
## 11 TRUE TRUE
## direct
## 1 Positive LFC
## 2 Positive LFC
## 3 Positive LFC
## 4 Negative LFC
## 5 Negative LFC
## 6 Negative LFC
## 7 Negative LFC
## 8 Negative LFC
## 9 Negative LFC
## 10 Negative LFC
## 11 Negative LFC# Make bar plot
dif.FLU.plot <- dif.FLU.df %>%
ggplot(aes(x = taxon, y = lfc_Former_landuseIndustrial, fill = direct)) +
geom_bar(stat = "identity", width = 0.7, color = "black",
position = position_dodge(width = 0.4)) +
geom_errorbar(aes(ymin = lfc_Former_landuseIndustrial - se_Former_landuseIndustrial,
ymax = lfc_Former_landuseIndustrial + se_Former_landuseIndustrial),
width = 0.2, position = position_dodge(0.05),
color = "black") +
labs(x = "Genus", y = "Log fold change") +
ggtitle(label = "Differentially abundant taxa",
subtitle="Prior land use: Industrial vs Agricultural") +
scale_fill_discrete(name = NULL) +
scale_color_discrete(name = NULL) +
theme_bw() +
theme(panel.grid.minor.y = element_blank(),
axis.text.x = element_text(size = 6, angle = 60, hjust = 1))
# Convert to plotly
dif.FLU.plot_plotly <- ggplotly(dif.FLU.plot, tooltip = "text")
dif.FLU.plot_plotly# Save the plot as an HTML file
htmlwidgets::saveWidget(dif.FLU.plot_plotly, "results/plotly/differential_abundant_taxa_FLU.html")# box plot function
diff_abund.box_plot <- function(taxon, env_factor, title) {
taxon_data <- psmelt(BacIndAgriGenera.rarefied) %>%
filter(Genus == taxon)
plot <- ggplot(taxon_data, aes_string(x = env_factor, y = "Abundance", fill = env_factor)) +
geom_violin() +
geom_jitter(width = 0.2, alpha = 0.5) +
labs(title = title, x = env_factor, y = "Abundance") +
theme_minimal() +
theme(legend.position = "none")
return(plot)
}# Extract the first positive and negative taxa
first_positive_taxa.FLU <- dif.FLU.df %>%
filter(direct == "Positive LFC") %>%
slice(1) %>%
pull(taxon)
first_negative_taxa.FLU <- dif.FLU.df %>%
filter(direct == "Negative LFC") %>%
slice(1) %>%
pull(taxon)
# Generate and plot the box plots
positive_plot.FLU <- diff_abund.box_plot(first_positive_taxa.FLU, "Former_landuse",
paste("Abundance of", first_positive_taxa.FLU, "by Former Land Use"))
negative_plot.FLU <- diff_abund.box_plot(first_negative_taxa.FLU, "Former_landuse",
paste("Abundance of", first_negative_taxa.FLU, "by Former Land Use"))
print(positive_plot.FLU)
print(negative_plot.FLU)
6.2 Region
# 2. Region
ancom.region = ancombc2(data = BacIndAgriGenera.rarefied2.merged_replica, tax_level = "Genus",
fix_formula = "Region",
group = "Region",
p_adj_method = "bonferroni", alpha = 0.05, pairwise = TRUE,
prv_cut = 0.10, lib_cut = 0,
n_cl = 4, verbose = TRUE)# Extract the results table
ancom.region.df <- ancom.region$res# Select only columns that we need
ancom.region.df <- ancom.region.df %>%
select(taxon, contains("Region"))# Select differentially abundant taxa
dif.region.df <- ancom.region.df %>%
filter(diff_RegionScotland & passed_ss_RegionScotland) %>%
arrange(desc(lfc_RegionScotland))sum(dif.region.df$lfc_RegionScotland > 1)## [1] 17# Filter out values that are not 2 fold change in both positive and negative direction
dif.region.df <- dif.region.df %>%
filter(lfc_RegionScotland >= 1 | lfc_RegionScotland <= -1)
# Create the 'direct' column to categorize LFC
dif.region.df <- dif.region.df %>%
mutate(
direct = ifelse(lfc_RegionScotland >= 1, "Positive LFC", "Negative LFC")
)
# Ensure direct factorized
dif.region.df$direct =
factor(dif.region.df$direct, levels = c("Positive LFC", "Negative LFC"))
dif.region.df$taxon =
factor(dif.region.df$taxon, levels = dif.region.df$taxon)dif.region.df## taxon lfc_RegionScotland se_RegionScotland W_RegionScotland
## 1 VBOC01 1.642275 0.1719702 9.549768
## 2 Palsa-189 1.511405 0.2475229 6.106121
## 3 Nemorincola 1.509296 0.1429288 10.559774
## 4 Flaviaesturariibacter 1.409598 0.2144276 6.573773
## 5 Pelomicrobium 1.391422 0.2058517 6.759339
## 6 2013-40CM-41-45 1.381802 0.1425396 9.694161
## 7 SCGC-AG-212-J23 1.264733 0.1737948 7.277162
## 8 Acidoferrum 1.208742 0.2791036 4.330801
## 9 Gp1-AA122 1.195852 0.2403979 4.974470
## 10 CADDZI01 1.182387 0.1965815 6.014746
## 11 UBA4416 1.178866 0.1712315 6.884633
## 12 Hypericibacter 1.158846 0.1634360 7.090519
## 13 Chitinophaga 1.156916 0.1527348 7.574666
## 14 Gp1-AA133 1.152696 0.2020446 5.705155
## 15 Palsa-187 1.114763 0.1918483 5.810648
## 16 RBG-16-71-46 1.047354 0.2228241 4.700364
## 17 UBA923 1.022077 0.1429087 7.151960
## 18 Actinoplanes -1.045388 0.1957841 -5.339495
## 19 VFJQ01 -1.048863 0.2197708 -4.772531
## 20 Methyloligella -1.096814 0.1415061 -7.751000
## 21 Luteitalea -1.129800 0.2329756 -4.849433
## 22 SHVA01 -1.143939 0.2414934 -4.736937
## 23 HRBIN40 -1.168912 0.2587924 -4.516792
## 24 Streptomyces_G_399870 -1.194410 0.1707691 -6.994297
## 25 Povalibacter -1.249674 0.2262127 -5.524331
## 26 Parviterribacter -1.351487 0.1549108 -8.724291
## 27 Microlunatus_A_391020 -1.353893 0.2000433 -6.768001
## 28 Solirubrobacter -1.405081 0.2846144 -4.936789
## 29 AC-51 -1.465112 0.2722885 -5.380735
## 30 VFJN01 -1.559506 0.2350581 -6.634555
## 31 Nocardioides_A_392796 -1.610780 0.2123895 -7.584086
## 32 Kribbella -1.665643 0.2004193 -8.310791
## 33 JAAAPF01 -1.705734 0.1759461 -9.694642
## 34 GMQP-bins7 -1.734388 0.2663932 -6.510632
## 35 Ohtaekwangia -1.861354 0.2649990 -7.024001
## p_RegionScotland q_RegionScotland diff_RegionScotland
## 1 1.805978e-08 7.476749e-06 TRUE
## 2 1.605755e-07 6.647824e-05 TRUE
## 3 1.678942e-10 6.950821e-08 TRUE
## 4 1.197430e-07 4.957361e-05 TRUE
## 5 8.640151e-09 3.577023e-06 TRUE
## 6 6.912439e-05 2.861750e-02 TRUE
## 7 6.712846e-09 2.779118e-06 TRUE
## 8 6.489750e-05 2.686757e-02 TRUE
## 9 7.015954e-06 2.904605e-03 TRUE
## 10 1.093486e-05 4.527032e-03 TRUE
## 11 2.134447e-07 8.836609e-05 TRUE
## 12 3.681885e-06 1.524300e-03 TRUE
## 13 4.808854e-10 1.990865e-07 TRUE
## 14 6.263632e-07 2.593143e-04 TRUE
## 15 3.424948e-07 1.417928e-04 TRUE
## 16 2.298778e-05 9.516941e-03 TRUE
## 17 4.919953e-06 2.036860e-03 TRUE
## 18 4.921712e-06 2.037589e-03 TRUE
## 19 1.741330e-05 7.209107e-03 TRUE
## 20 1.114786e-04 4.615212e-02 TRUE
## 21 2.690206e-05 1.113745e-02 TRUE
## 22 1.663815e-05 6.888195e-03 TRUE
## 23 3.542722e-05 1.466687e-02 TRUE
## 24 2.996486e-09 1.240545e-06 TRUE
## 25 1.324347e-06 5.482796e-04 TRUE
## 26 1.100480e-05 4.555987e-03 TRUE
## 27 5.075079e-08 2.101083e-05 TRUE
## 28 8.904792e-06 3.686584e-03 TRUE
## 29 2.884837e-06 1.194323e-03 TRUE
## 30 4.362516e-08 1.806081e-05 TRUE
## 31 3.056924e-10 1.265567e-07 TRUE
## 32 4.121611e-11 1.706347e-08 TRUE
## 33 2.623744e-05 1.086230e-02 TRUE
## 34 1.935353e-08 8.012362e-06 TRUE
## 35 4.987584e-09 2.064860e-06 TRUE
## passed_ss_RegionScotland direct
## 1 TRUE Positive LFC
## 2 TRUE Positive LFC
## 3 TRUE Positive LFC
## 4 TRUE Positive LFC
## 5 TRUE Positive LFC
## 6 TRUE Positive LFC
## 7 TRUE Positive LFC
## 8 TRUE Positive LFC
## 9 TRUE Positive LFC
## 10 TRUE Positive LFC
## 11 TRUE Positive LFC
## 12 TRUE Positive LFC
## 13 TRUE Positive LFC
## 14 TRUE Positive LFC
## 15 TRUE Positive LFC
## 16 TRUE Positive LFC
## 17 TRUE Positive LFC
## 18 TRUE Negative LFC
## 19 TRUE Negative LFC
## 20 TRUE Negative LFC
## 21 TRUE Negative LFC
## 22 TRUE Negative LFC
## 23 TRUE Negative LFC
## 24 TRUE Negative LFC
## 25 TRUE Negative LFC
## 26 TRUE Negative LFC
## 27 TRUE Negative LFC
## 28 TRUE Negative LFC
## 29 TRUE Negative LFC
## 30 TRUE Negative LFC
## 31 TRUE Negative LFC
## 32 TRUE Negative LFC
## 33 TRUE Negative LFC
## 34 TRUE Negative LFC
## 35 TRUE Negative LFC# Make bar plot
dif.region.plot <- dif.region.df %>%
ggplot(aes(x = taxon, y = lfc_RegionScotland, fill = direct)) +
geom_bar(stat = "identity", width = 0.7, color = "black",
position = position_dodge(width = 0.4)) +
geom_errorbar(aes(ymin = lfc_RegionScotland - se_RegionScotland,
ymax = lfc_RegionScotland + se_RegionScotland),
width = 0.2, position = position_dodge(0.05),
color = "black") +
labs(x = "Genus", y = "Log fold change") +
ggtitle(label = "Differentially abundant taxa",
subtitle="Region: Scotland vs England") +
scale_fill_discrete(name = NULL) +
scale_color_discrete(name = NULL) +
theme_bw() +
theme(panel.grid.minor.y = element_blank(),
axis.text.x = element_text(size = 6, angle = 60, hjust = 1))
# Convert to plotly
dif.region.plot_plotly <- ggplotly(dif.region.plot, tooltip = "text")
dif.region.plot_plotly# Save the plot as an HTML file
htmlwidgets::saveWidget(dif.region.plot_plotly, "results/plotly/differential_abundant_taxa_Region.html")# Extract the first positive and negative taxa
first_positive_taxa.region <- dif.region.df %>%
filter(direct == "Positive LFC") %>%
slice(1) %>%
pull(taxon)
first_negative_taxa.region <- dif.region.df %>%
filter(direct == "Negative LFC") %>%
slice(1) %>%
pull(taxon)
# Generate and plot the box plots
positive_plot.region <- diff_abund.box_plot(first_positive_taxa.region, "Region",
paste("Abundance of", first_positive_taxa.region, "by Region"))
negative_plot.region <- diff_abund.box_plot(first_negative_taxa.region, "Region",
paste("Abundance of", first_negative_taxa.region, "by Region"))
print(positive_plot.region)
print(negative_plot.region)
7 Ordination
7.1 NMDS
# Compute unweighted UniFrac distance
unifrac_dist <- UniFrac(BacIndAgriGenera.rarefied2.merged_replica, weighted = FALSE)
# Perform NMDS ordination
ordination.NMDS <- ordinate(BacIndAgriGenera.rarefied2.merged_replica, method = "NMDS", distance = unifrac_dist)
## For 3d plotting
# NMDS ordination
ordination.NMDS.3d <- metaMDS(unifrac_dist, k=3)
# Extract NMDS samples scores (the position of each sample in the ordination space)
nmds_scores <- scores(ordination.NMDS.3d, display="sites")
# Convert to data frame and add sample metadata
nmds_data <- as.data.frame(nmds_scores)###Region
# 2d plot
ordination_NMDS.Region <- plot_ordination(BacIndAgriGenera.rarefied2.merged_replica, ordination.NMDS, type="samples", color="Region") +
ggtitle("samples ordination:",
subtitle="NMDS with unweighted UniFrac distances")+
geom_point(size=3) + stat_ellipse() +
theme_bw() +
theme(legend.position = "bottom") +
guides(color=guide_legend(ncol=3))
print(ordination_NMDS.Region)
nmds_data$Region <- sample_data(BacIndAgriGenera.rarefied2.merged_replica)$Region
# 3D plot
ordination_NMDS.Region.3d <- plot_ly(nmds_data, x = ~NMDS1, y = ~NMDS2, z = ~NMDS3,
color = ~Region, colors = c("red", "blue"),
type = 'scatter3d', mode = 'markers',
marker = list(size = 5))
ordination_NMDS.Region.3d <- ordination_NMDS.Region.3d %>% layout(title = "NMDS with Unweighted UniFrac Distances",
scene = list(xaxis = list(title = 'NMDS1'),
yaxis = list(title = 'NMDS2'),
zaxis = list(title = 'NMDS3')))
ordination_NMDS.Region.3d###FLU
ordination_NMDS.FLU <- plot_ordination(BacIndAgriGenera.rarefied2.merged_replica, ordination.NMDS, type="samples", color="Former_landuse") +
ggtitle("samples ordination:",
subtitle="NMDS with unweighted UniFrac distances")+
geom_point(size=3) + stat_ellipse() +
theme_bw() +
theme(legend.position = "bottom") +
guides(color=guide_legend(ncol=3))
print(ordination_NMDS.FLU)
nmds_data$Former_landuse <- sample_data(BacIndAgriGenera.rarefied2.merged_replica)$Former_landuse
# 3D plot
ordination_NMDS.FLU.3d <- plot_ly(nmds_data, x = ~NMDS1, y = ~NMDS2, z = ~NMDS3,
color = ~Former_landuse, colors = c("red", "blue"),
type = 'scatter3d', mode = 'markers',
marker = list(size = 5))
ordination_NMDS.FLU.3d <- ordination_NMDS.FLU.3d %>% layout(title = "NMDS with Unweighted UniFrac Distances",
scene = list(xaxis = list(title = 'NMDS1'),
yaxis = list(title = 'NMDS2'),
zaxis = list(title = 'NMDS3')))
ordination_NMDS.FLU.3d7.2 RDA
# Convert BacIndAgriGenera.rarefied2.merged_replica to data frame for abundance and environmental data
abund_data <- as.data.frame(t(otu_table(BacIndAgriGenera.rarefied2.merged_replica)))
env_data <- data.frame(as(sample_data(BacIndAgriGenera.rarefied2.merged_replica), "data.frame"))
# Select variables
env_data <- env_data[, c("Region", "Former_landuse", "Woodland_age", "pH", "Conductivity")]
# Make sure pH, age, and conductivity are numerical
env_data$Woodland_age <- as.numeric(as.character(env_data$Woodland_age))
env_data$Conductivity <- as.numeric(as.character(env_data$Conductivity))
env_data$pH <- as.numeric(as.character(env_data$pH))
# Remove NA values
env_data <- na.omit(env_data)
# Factorize categorical columns
env_data$Region <- as.factor(env_data$Region)
env_data$Former_landuse <- as.factor(env_data$Former_landuse)
# Select row names in abund_data that match env_data row names
abund_data <- abund_data[rownames(abund_data) %in% rownames(env_data), ]
# Perform RDA, constraining all env_data variables
ordination.RDA <- rda(abund_data ~ pH + Conductivity, data = env_data)RDA.R2adj <- RsquareAdj(ordination.RDA)$adj.r.squared
RDA.R2adj## [1] 0.356577# # Enable the r-universe repo
# options(repos = c(
# fawda123 = 'https://fawda123.r-universe.dev',
# CRAN = 'https://cloud.r-project.org'))
#
# # Install ggord
# install.packages('ggord')library(ggord)ggord(ordination.RDA, env_data$Region, ext=0.8, txt=4, grp_title="Region", xlims=c(-1,1), ylims=c(-1,1), alpha_el=0.6, size=1.5, poly=FALSE,
vectyp="longdash", obslab=FALSE, ptslab=FALSE) 
ggord(ordination.RDA, env_data$Former_landuse, ext=0.8, txt=4, grp_title="Former_landuse", xlims=c(-1,1), ylims=c(-1,1), poly=FALSE, size=1.5,
vectyp="longdash", obslab=FALSE, ptslab=FALSE) 
7.3 cCA
# Perform CCA
ordination.cCA <- cca(abund_data ~ pH + Conductivity, data = env_data)cCA.R2adj <- RsquareAdj(ordination.cCA)$adj.r.squared
cCA.R2adj## [1] 0.2424402ggord(ordination.cCA, env_data$Region, ext=0.8, txt=4, grp_title="Region", xlims=c(-1,1), ylims=c(-1,1), alpha_el=0.6, size=1.5,
vectyp="longdash", obslab=FALSE) 
ggord(ordination.cCA, env_data$Former_landuse, ext=0.8, txt=4, grp_title="Former_landuse", xlims=c(-1,1), ylims=c(-1,1), alpha_el=0.6, size=1.5,
vectyp="longdash", obslab=FALSE) 
8 Beta Diversity (nested design)
# Re transpose abund_data
otu_table_matrix <- as.matrix(t(abund_data))
# Create a new phyloseq object with the filtered samples (based on env_data, no NaNs)
BacIndAgriGenera.rarefied2.merged.filtered <- phyloseq(
otu_table(otu_table_matrix, taxa_are_rows = TRUE),
sample_data(env_data),
tax_table(BacIndAgriGenera.rarefied2.merged_replica),
phy_tree(BacIndAgriGenera.rarefied2.merged_replica)
)# Calculate beta diversity metrics
weighted_unifrac_dist.filtered <- UniFrac(BacIndAgriGenera.rarefied2.merged.filtered, weighted = TRUE)
unweighted_unifrac_dist.filtered <- UniFrac(BacIndAgriGenera.rarefied2.merged.filtered, weighted = FALSE)
bray_dist.filtered <- distance(BacIndAgriGenera.rarefied2.merged.filtered, method="bray")adjus_p_value <- function(adonis_result){
# Extract p-values
p_values <- adonis_result$`Pr(>F)`
# Adjust the p-values
p_adjusted <- p.adjust(p_values)
# Add the adjusted p-values to adonis2 table
adonis_result$`Pr(>F)_adjusted` <- c(p_adjusted)
return(adonis_result)
}library(scales)
nested_model.pie <- function(nested_model, title){
# Extract the explained variance (R2)
explained_variance <- nested_model$R2
# Remove the NA values associated with the rest of table info
p_values <- nested_model$`Pr(>F)`
explained_variance <- explained_variance[!is.na(p_values)]
# Names of the factors
factors <- rownames(nested_model)[1:length(explained_variance)]
# Create a data frame for plotting
plot_data <- data.frame(
Factor = factors,
ExplainedVariance = explained_variance
)
# Ensure the levels of Factor are in the same order as the explained variance
plot_data$Factor <- factor(plot_data$Factor, levels = plot_data$Factor)
# Custom legends
plot_data$LegendLabel <- paste0(plot_data$Factor,
" (", percent(plot_data$ExplainedVariance), ")")
# Create the pie chart
pie_chart <- ggplot(plot_data, aes(x = "", y = ExplainedVariance, fill = Factor)) +
geom_bar(width = 1, stat = "identity") +
coord_polar(theta = "y") +
scale_y_continuous(labels = percent_format()) +
theme_void() +
theme(legend.position = "right") +
labs(title = title, fill = "Factor") +
scale_fill_manual(values = scales::brewer_pal(palette = "Set3")(length(plot_data$Factor)),
labels = plot_data$LegendLabel)
print(pie_chart)
}pH / conductivity / FLU / woodland age / region
8.1 Weighted UniFrac
nested_model.weighted_unifrac <- adonis2(weighted_unifrac_dist.filtered ~ Region/Former_landuse
+ pH + Woodland_age,
data = env_data, strata = env_data$Region, permutations=1000)
nested_model.weighted_unifrac## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = weighted_unifrac_dist.filtered ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
## Df SumOfSqs R2 F Pr(>F)
## Region 1 0.85916 0.31705 46.3139 0.000999 ***
## pH 1 0.58442 0.21567 31.5038 0.000999 ***
## Woodland_age 1 0.03752 0.01385 2.0225 0.088911 .
## Region:Former_landuse 2 0.24555 0.09061 6.6182 0.000999 ***
## Residual 53 0.98319 0.36282
## Total 58 2.70984 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adjus_p_value(nested_model.weighted_unifrac)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = weighted_unifrac_dist.filtered ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
## Df SumOfSqs R2 F Pr(>F) Pr(>F)_adjusted
## Region 1 0.85916 0.31705 46.3139 0.000999 0.003996 **
## pH 1 0.58442 0.21567 31.5038 0.000999 0.003996 **
## Woodland_age 1 0.03752 0.01385 2.0225 0.088911 0.088911 .
## Region:Former_landuse 2 0.24555 0.09061 6.6182 0.000999 0.003996 **
## Residual 53 0.98319 0.36282
## Total 58 2.70984 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nested_model.pie(nested_model.weighted_unifrac, "Explained Variance of weighted UniFrac data by environmental factors")
8.2 Unweighted UniFrac
nested_model.unweighted_unifrac <- adonis2(unweighted_unifrac_dist.filtered ~ Region/Former_landuse
+ pH + Woodland_age,
data = env_data, strata = env_data$Region, permutations=1000)
nested_model.unweighted_unifrac## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = unweighted_unifrac_dist.filtered ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
## Df SumOfSqs R2 F Pr(>F)
## Region 1 0.7242 0.19127 17.4608 0.000999 ***
## pH 1 0.4694 0.12396 11.3164 0.000999 ***
## Woodland_age 1 0.0719 0.01900 1.7345 0.049950 *
## Region:Former_landuse 2 0.3226 0.08520 3.8889 0.000999 ***
## Residual 53 2.1982 0.58057
## Total 58 3.7863 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adjus_p_value(nested_model.unweighted_unifrac)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = unweighted_unifrac_dist.filtered ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
## Df SumOfSqs R2 F Pr(>F) Pr(>F)_adjusted
## Region 1 0.7242 0.19127 17.4608 0.000999 0.003996 **
## pH 1 0.4694 0.12396 11.3164 0.000999 0.003996 **
## Woodland_age 1 0.0719 0.01900 1.7345 0.049950 0.049950 *
## Region:Former_landuse 2 0.3226 0.08520 3.8889 0.000999 0.003996 **
## Residual 53 2.1982 0.58057
## Total 58 3.7863 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nested_model.pie(nested_model.unweighted_unifrac, "Explained Variance of unweighted UniFrac data by environmental factors")
8.3 Bray curtis
nested_model.bray <- adonis2(bray_dist.filtered ~ Region/Former_landuse
+ pH + Woodland_age,
data = env_data, strata = env_data$Region, permutations=1000)
nested_model.bray## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = bray_dist.filtered ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
## Df SumOfSqs R2 F Pr(>F)
## Region 1 1.1178 0.22339 28.1794 0.000999 ***
## pH 1 1.2356 0.24694 31.1505 0.000999 ***
## Woodland_age 1 0.1036 0.02070 2.6115 0.027972 *
## Region:Former_landuse 2 0.4445 0.08883 5.6028 0.000999 ***
## Residual 53 2.1023 0.42014
## Total 58 5.0038 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adjus_p_value(nested_model.bray)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = bray_dist.filtered ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
## Df SumOfSqs R2 F Pr(>F) Pr(>F)_adjusted
## Region 1 1.1178 0.22339 28.1794 0.000999 0.003996 **
## pH 1 1.2356 0.24694 31.1505 0.000999 0.003996 **
## Woodland_age 1 0.1036 0.02070 2.6115 0.027972 0.027972 *
## Region:Former_landuse 2 0.4445 0.08883 5.6028 0.000999 0.003996 **
## Residual 53 2.1023 0.42014
## Total 58 5.0038 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nested_model.pie(nested_model.bray, "Explained Variance of bray curtis by environmental factors")
8.4 Abundance data
nested_model.abund <- adonis2(abund_data ~ Region/Former_landuse
+ pH + Woodland_age,
data = env_data, strata = env_data$Region, permutations=1000)
nested_model.abund## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = abund_data ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
## Df SumOfSqs R2 F Pr(>F)
## Region 1 1.1178 0.22339 28.1794 0.000999 ***
## pH 1 1.2356 0.24694 31.1505 0.000999 ***
## Woodland_age 1 0.1036 0.02070 2.6115 0.032967 *
## Region:Former_landuse 2 0.4445 0.08883 5.6028 0.000999 ***
## Residual 53 2.1023 0.42014
## Total 58 5.0038 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adjus_p_value(nested_model.abund)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = abund_data ~ Region/Former_landuse + pH + Woodland_age, data = env_data, permutations = 1000, strata = env_data$Region)
## Df SumOfSqs R2 F Pr(>F) Pr(>F)_adjusted
## Region 1 1.1178 0.22339 28.1794 0.000999 0.003996 **
## pH 1 1.2356 0.24694 31.1505 0.000999 0.003996 **
## Woodland_age 1 0.1036 0.02070 2.6115 0.032967 0.032967 *
## Region:Former_landuse 2 0.4445 0.08883 5.6028 0.000999 0.003996 **
## Residual 53 2.1023 0.42014
## Total 58 5.0038 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nested_model.pie(nested_model.abund, "Explained Variance of abundance data by environmental factors")